This tutorial is part of viskex.

import typing

import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import ufl

import viskex

Generate and plot meshes

interval = dolfinx.mesh.create_unit_interval(mpi4py.MPI.COMM_WORLD, 6)
square_tria = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 3, 3, dolfinx.mesh.CellType.triangle)
square_quad = dolfinx.mesh.create_unit_square(mpi4py.MPI.COMM_WORLD, 3, 3, dolfinx.mesh.CellType.quadrilateral)
cube_tetra = dolfinx.mesh.create_unit_cube(mpi4py.MPI.COMM_WORLD, 3, 3, 3, dolfinx.mesh.CellType.tetrahedron)
cube_hexa = dolfinx.mesh.create_unit_cube(mpi4py.MPI.COMM_WORLD, 3, 3, 3, dolfinx.mesh.CellType.hexahedron)

viskex.dolfinx.plot_mesh(interval)

viskex.dolfinx.plot_mesh(interval, dim=1)

viskex.dolfinx.plot_mesh(interval, dim=0)

viskex.dolfinx.plot_mesh(square_tria)

viskex.dolfinx.plot_mesh(square_tria, dim=2)

viskex.dolfinx.plot_mesh(square_tria, dim=1)

viskex.dolfinx.plot_mesh(square_tria, dim=0)

viskex.dolfinx.plot_mesh(square_quad)

viskex.dolfinx.plot_mesh(square_quad, dim=2)

viskex.dolfinx.plot_mesh(square_quad, dim=1)

viskex.dolfinx.plot_mesh(square_quad, dim=0)

viskex.dolfinx.plot_mesh(cube_tetra)

viskex.dolfinx.plot_mesh(cube_tetra, dim=3)

viskex.dolfinx.plot_mesh(cube_tetra, dim=2)

viskex.dolfinx.plot_mesh(cube_tetra, dim=1)

viskex.dolfinx.plot_mesh(cube_tetra, dim=0)

viskex.dolfinx.plot_mesh(cube_hexa)

viskex.dolfinx.plot_mesh(cube_hexa, dim=3)

viskex.dolfinx.plot_mesh(cube_hexa, dim=2)

viskex.dolfinx.plot_mesh(cube_hexa, dim=1)

viskex.dolfinx.plot_mesh(cube_hexa, dim=0)

Generate and plot subdomains

def mark_subdomains(mesh: dolfinx.mesh.Mesh) -> dolfinx.mesh.MeshTags:
    """Mark left and right subdomains in a given mesh with values 1 and 2, respectively."""
    def left_subdomain(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
        """Condition that defines the left subdomain."""
        return x[0] <= 1.0 / 3.0  # type: ignore[no-any-return]

    def right_subdomain(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
        """Condition that defines the right subdomain."""
        return x[0] >= 2.0 / 3.0  # type: ignore[no-any-return]

    subdomains_entities = dict()
    subdomains_values = dict()
    for (subdomain, subdomain_id) in zip((left_subdomain, right_subdomain), (1, 2)):
        subdomains_entities[subdomain_id] = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, subdomain)
        subdomains_values[subdomain_id] = np.full(subdomains_entities[subdomain_id].shape, subdomain_id, dtype=np.int32)

    subdomains_entities_unsorted = np.hstack(list(subdomains_entities.values()))
    subdomains_values_unsorted = np.hstack(list(subdomains_values.values()))
    subdomains_entities_argsort = np.argsort(subdomains_entities_unsorted)
    subdomains_entities_sorted = subdomains_entities_unsorted[subdomains_entities_argsort]
    subdomains_values_sorted = subdomains_values_unsorted[subdomains_entities_argsort]
    subdomains = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, subdomains_entities_sorted, subdomains_values_sorted)
    return subdomains

interval_subdomains = mark_subdomains(interval)
square_tria_subdomains = mark_subdomains(square_tria)
square_quad_subdomains = mark_subdomains(square_quad)
cube_tetra_subdomains = mark_subdomains(cube_tetra)
cube_hexa_subdomains = mark_subdomains(cube_hexa)

viskex.dolfinx.plot_mesh_entities(
    interval, 1, "subdomains", interval_subdomains.indices[interval_subdomains.values == 2])

viskex.dolfinx.plot_mesh_entities(
    interval, 1, "subdomains", interval_subdomains.indices, interval_subdomains.values)

viskex.dolfinx.plot_mesh_tags(interval, interval_subdomains, "subdomains")

viskex.dolfinx.plot_mesh_entities(
    square_tria, 2, "subdomains", square_tria_subdomains.indices[square_tria_subdomains.values == 2])

viskex.dolfinx.plot_mesh_entities(
    square_tria, 2, "subdomains", square_tria_subdomains.indices, square_tria_subdomains.values)

viskex.dolfinx.plot_mesh_tags(square_tria, square_tria_subdomains, "subdomains")

viskex.dolfinx.plot_mesh_tags(square_quad, square_quad_subdomains, "subdomains")

viskex.dolfinx.plot_mesh_tags(cube_tetra, cube_tetra_subdomains, "subdomains")

viskex.dolfinx.plot_mesh_tags(cube_hexa, cube_hexa_subdomains, "subdomains")

Generate and plot boundaries

def mark_boundaries(mesh: dolfinx.mesh.Mesh, subdomains: dolfinx.mesh.MeshTags) -> dolfinx.mesh.MeshTags:
    """
    Mark internal and boundary facets in a given mesh with four different values.

    Internal facets of left and right subdomains are associated with values 1 and 2, respectively.
    Furthermore, boundary facets on the left and right boundaries are associated with values 3 and 4,
    respectively.
    """
    def left(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
        """Condition that defines the left boundary."""
        return abs(x[0] - 0.) < np.finfo(float).eps  # type: ignore[no-any-return]

    def right(x: np.typing.NDArray[np.float64]) -> np.typing.NDArray[np.bool_]:
        """Condition that defines the right boundary."""
        return abs(x[0] - 1.) < np.finfo(float).eps  # type: ignore[no-any-return]

    boundaries_entities = dict()
    boundaries_values = dict()
    for (boundary, boundary_id) in zip((left, right), (3, 4)):
        boundaries_entities[boundary_id] = dolfinx.mesh.locate_entities_boundary(
            mesh, mesh.topology.dim - 1, boundary)
        boundaries_values[boundary_id] = np.full(
            boundaries_entities[boundary_id].shape, boundary_id, dtype=np.int32)

    boundaries_entities_unsorted = np.hstack(list(boundaries_entities.values()))
    boundaries_values_unsorted = np.hstack(list(boundaries_values.values()))
    boundaries_entities_argsort = np.argsort(boundaries_entities_unsorted)
    boundaries_entities_sorted = boundaries_entities_unsorted[boundaries_entities_argsort]
    boundaries_values_sorted = boundaries_values_unsorted[boundaries_entities_argsort]
    boundaries = dolfinx.mesh.meshtags(
        mesh, mesh.topology.dim - 1,
        boundaries_entities_sorted, boundaries_values_sorted)

    cell_to_facets_connectivity = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim - 1)
    new_facets_indices = list()
    new_facets_values = list()
    for cell in range(cell_to_facets_connectivity.num_nodes):
        facets = cell_to_facets_connectivity.links(cell)
        subdomain_index = np.where(subdomains.indices == cell)[0]
        if subdomain_index.size > 0:
            assert subdomain_index.size == 1
            subdomain_value = subdomains.values[subdomain_index[0]]
            if subdomain_value in (1, 2):
                for facet in facets:
                    if facet not in boundaries.indices and facet not in new_facets_indices:
                        new_facets_indices.append(facet)
                        new_facets_values.append(subdomain_value)

    boundaries_and_interfaces_entities_unsorted = np.hstack((boundaries.indices, new_facets_indices)).astype(np.int32)
    boundaries_and_interfaces_values_unsorted = np.hstack((boundaries.values, new_facets_values)).astype(np.int32)
    boundaries_and_interfaces_entities_argsort = np.argsort(boundaries_and_interfaces_entities_unsorted)
    boundaries_and_interfaces_entities_sorted = boundaries_and_interfaces_entities_unsorted[
        boundaries_and_interfaces_entities_argsort]
    boundaries_and_interfaces_values_sorted = boundaries_and_interfaces_values_unsorted[
        boundaries_and_interfaces_entities_argsort]
    boundaries_and_interfaces = dolfinx.mesh.meshtags(
        mesh, mesh.topology.dim - 1,
        boundaries_and_interfaces_entities_sorted, boundaries_and_interfaces_values_sorted)

    return boundaries_and_interfaces

interval_boundaries = mark_boundaries(interval, interval_subdomains)
square_tria_boundaries = mark_boundaries(square_tria, square_tria_subdomains)
square_quad_boundaries = mark_boundaries(square_quad, square_quad_subdomains)
cube_tetra_boundaries = mark_boundaries(cube_tetra, cube_tetra_subdomains)
cube_hexa_boundaries = mark_boundaries(cube_hexa, cube_hexa_subdomains)

viskex.dolfinx.plot_mesh_entities(
    interval, 0, "boundaries", interval_boundaries.indices[interval_boundaries.values == 2])

viskex.dolfinx.plot_mesh_entities(
    interval, 0, "boundaries", interval_boundaries.indices, interval_boundaries.values)

viskex.dolfinx.plot_mesh_tags(interval, interval_boundaries, "boundaries")

viskex.dolfinx.plot_mesh_entities(
    square_tria, 1, "boundaries", square_tria_boundaries.indices[square_tria_boundaries.values == 2])

viskex.dolfinx.plot_mesh_entities(
    square_tria, 1, "boundaries", square_tria_boundaries.indices, square_tria_boundaries.values)

viskex.dolfinx.plot_mesh_tags(square_tria, square_tria_boundaries, "boundaries")

viskex.dolfinx.plot_mesh_tags(square_quad, square_quad_boundaries, "boundaries")

viskex.dolfinx.plot_mesh_tags(cube_tetra, cube_tetra_boundaries, "boundaries")

viskex.dolfinx.plot_mesh_tags(cube_hexa, cube_hexa_boundaries, "boundaries")

Interpolate and plot scalar functions

def prepare_scalar_field_cases(  # type: ignore[no-any-unimported]
    mesh: dolfinx.mesh.Mesh,
    expression: typing.Callable[
        [typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]],
        typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]
    ]
) -> typing.Tuple[dolfinx.fem.Function, typing.Tuple[ufl.core.expr.Expr, dolfinx.fem.FunctionSpace]]:
    """Prepare scalar field cases."""
    scalar_function_space = dolfinx.fem.FunctionSpace(mesh, ("CG", 2))
    scalar_field = dolfinx.fem.Function(scalar_function_space)
    scalar_field.interpolate(expression)
    scalar_field_ufl = expression(ufl.SpatialCoordinate(mesh))
    return scalar_field, (scalar_field_ufl, scalar_function_space)

interval_scalar_field, interval_scalar_field_ufl = prepare_scalar_field_cases(
    interval, lambda x: x[0]**3)
square_tria_scalar_field, square_tria_scalar_field_ufl = prepare_scalar_field_cases(
    square_tria, lambda x: x[0]**3 + x[1]**2)
square_quad_scalar_field, square_quad_scalar_field_ufl = prepare_scalar_field_cases(
    square_quad, lambda x: x[0]**3 + x[1]**2)
cube_tetra_scalar_field, cube_tetra_scalar_field_ufl = prepare_scalar_field_cases(
    cube_tetra, lambda x: x[0]**3 + x[1]**2 + x[2]**4)
cube_hexa_scalar_field, cube_hexa_scalar_field_ufl = prepare_scalar_field_cases(
    cube_hexa, lambda x: x[0]**3 + x[1]**2 + x[2]**4)

viskex.dolfinx.plot_scalar_field(interval_scalar_field, "scalar")

viskex.dolfinx.plot_scalar_field(interval_scalar_field_ufl, "scalar")

viskex.dolfinx.plot_scalar_field(square_tria_scalar_field, "scalar")

viskex.dolfinx.plot_scalar_field(square_tria_scalar_field, "scalar", warp_factor=1.0)

viskex.dolfinx.plot_scalar_field(square_tria_scalar_field_ufl, "scalar")

viskex.dolfinx.plot_scalar_field(square_quad_scalar_field, "scalar")

viskex.dolfinx.plot_scalar_field(square_quad_scalar_field, "scalar", warp_factor=1.0)

viskex.dolfinx.plot_scalar_field(square_quad_scalar_field_ufl, "scalar")

viskex.dolfinx.plot_scalar_field(cube_tetra_scalar_field, "scalar")

viskex.dolfinx.plot_scalar_field(cube_tetra_scalar_field_ufl, "scalar")

viskex.dolfinx.plot_scalar_field(cube_hexa_scalar_field, "scalar")

viskex.dolfinx.plot_scalar_field(cube_hexa_scalar_field_ufl, "scalar")

Interpolate and plot vector functions

def prepare_vector_field_cases(  # type: ignore[no-any-unimported]
    mesh: dolfinx.mesh.Mesh,
    expression: typing.Callable[
        [typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]],
        typing.Union[np.typing.NDArray[np.float64], ufl.core.expr.Expr]
    ]
) -> typing.Tuple[dolfinx.fem.Function, typing.Tuple[ufl.core.expr.Expr, dolfinx.fem.FunctionSpace]]:
    """Prepare vector field cases."""
    vector_function_space = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2))
    vector_field = dolfinx.fem.Function(vector_function_space)
    vector_field.interpolate(expression)
    vector_field_ufl = ufl.as_vector(expression(ufl.SpatialCoordinate(mesh)))
    return vector_field, (vector_field_ufl, vector_function_space)

square_tria_vector_field, square_tria_vector_field_ufl = prepare_vector_field_cases(
    square_tria, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
square_quad_vector_field, square_quad_vector_field_ufl = prepare_vector_field_cases(
    square_quad, lambda x: (x[0]**3 + x[1]**2, x[0]**5 + x[1]**4))
cube_tetra_vector_field, cube_tetra_vector_field_ufl = prepare_vector_field_cases(
    cube_tetra, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))
cube_hexa_vector_field, cube_hexa_vector_field_ufl = prepare_vector_field_cases(
    cube_hexa, lambda x: (x[0]**3 + x[1]**2 + x[2]**4, x[0]**6 + x[1]**5 + x[2]**7, x[0]**9 + x[1]**8 + x[2]**10))

viskex.dolfinx.plot_vector_field(square_tria_vector_field, "vector")

viskex.dolfinx.plot_vector_field(square_tria_vector_field, "vector", glyph_factor=0.1)

viskex.dolfinx.plot_vector_field(square_tria_vector_field, "vector", warp_factor=1.0)

viskex.dolfinx.plot_vector_field(square_tria_vector_field_ufl, "vector")

viskex.dolfinx.plot_vector_field(square_quad_vector_field, "vector")

viskex.dolfinx.plot_vector_field(square_quad_vector_field, "vector", glyph_factor=0.1)

viskex.dolfinx.plot_vector_field(square_quad_vector_field, "vector", warp_factor=1.0)

viskex.dolfinx.plot_vector_field(cube_tetra_vector_field, "vector")

viskex.dolfinx.plot_vector_field(cube_tetra_vector_field, "vector", glyph_factor=0.1)

viskex.dolfinx.plot_vector_field(cube_tetra_vector_field, "vector", warp_factor=1.0)

viskex.dolfinx.plot_vector_field(cube_hexa_vector_field, "vector")

viskex.dolfinx.plot_vector_field(cube_hexa_vector_field, "vector", glyph_factor=0.1)

viskex.dolfinx.plot_vector_field(cube_hexa_vector_field, "vector", warp_factor=1.0)